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ABSTRACT 


A  general-purpose  viscous  flow  code,  DTNS3D,  has  been  used  to  compute  the 
strength  and  trajectories  of  an  isolated  vortex  and  a  vortex  pair.  Analytic  solutions 
and  previous  computations  are  used  to  judge  the  accuracy  of  the  three-dimensional 
steady-state  results.  It  is  possible  to  track,  quite  accurately,  the  decay  and  translation 
of  these  vortex  structures.  A  byproduct  of  this  study  is  some  information  on  the 
limitations  of  the  unsteady  crossflow  analogy  and  on  the  effect  of  axial  velocity 
variations  on  vortex  trajectories. 

ADMINISTRATIVE  INFORMATION 

This  report  is  submitted  in  partial  fulfillment  of  Task  2,  Hydrodynamic  Wake  Signatures  in 
the  Maneuvering  and  Control  Project,  RB23H16,  in  the  Submarine  Technology  Block  Program 
(ND3A/PE0602323N)  for  fiscal  year  1994.  The  work  described  herein  was  sponsored  by  the 
Office  of  Naval  Research  (ONR  334)  and  performed  by  the  Carderock  Division,  Naval  Surface 
Warfare  Center  under  Work  Unit  5060-435. 

INTRODUCTION 

Vortices  shed  from  objects  moving  through  air  and  water  as  well  as  from  fixed  objects  in 
flowing  fluids  are  among  the  most  striking  and  significant  features  of  such  flows.  These  flow 
strucmres  have  important  implications  in  such  areas  as  aircraft  safety,  sediment  transport,  and 
marine  vehicle  maneuvering  and  control.  The  abundant  effects  of  these  vortices  are  among  the 
reasons  that  much  effort  has  been  put  into  the  development  of  computer  codes  for  the  numerical 
solution  of  equations  governing  viscous  fluid  motions.  The  successful  construction  of  such 
codes  will  allow  the  prediction  of  vortex  generation  and  movement  and  will  lead  to  a  design 
capability  for  suppression  or  other  modification  of  vortex  flow  features.  Calculations  are 
already  being  carried  out  regularly  in  many  industries  for  flows  around  three-dimensional 
bodies.  At  Carderock  Division,  Naval  Surface  Warfare  Center  (NSWCCD)  it  is  now  quite 
common  for  viscous  flow  computations  to  be  carried  out  in  parallel  with  experimental  efforts 
aimed  at  understanding  and  improving  the  flow  about  naval  vehicles.  These  computations 
involve  vortices  shed  from  appendages,  and  indications  are  that  the  computer  codes  being  used 
can  track  the  vortices  at  least  in  a  qualitatively  correct  manner.  However,  a  large  part  of  the 
ongoing  validation  of  these  codes  must  be  a  careful  assessment  of  their  ability  to  accurately 
predict  vorticity  production,  concentration,  diffusion,  and  convection.  Only  when  such  an 
assessment  is  complete  will  it  be  known  to  what  extent  and  how  these  codes  can  be  used  as 
reliable  research  and  design  tools. 

For  these  reasons,  an  effon  was  recently  initiated  at  NSWCCD  to  validate  the  capability  of 
existing  methods  and  computer  codes  to  accurately  predict  the  generation  and  movement  of 
vortices.  This  is  a  difficult  task.  Reliable  experimental  data  are  scarce.  Most  flows  of  interest 
are  mrbulent.  When  computational  results  are  compared  to  measured  data  it  is  usually 
impossible  to  quantitatively  divide  the  discrepancies  into  such  categories  as  numerical  error, 
measurement  error,  and  shortcomings  of  turbulence  models.  To  render  this  validation  process 
more  manageable,  the  effort  has  been  divided  into  the  nearfield  task  of  vorticity  production  and 
vortex  formation  and  the  farfield  task  of  vortex  tracking.  The  work  being  reported  falls  into  the 
farfield  category. 

This  report  is  one  of  a  series  describing  this  validation  effort.  In  previous  documentation, 
Ohring  and  Lugt'  showed  that  a  careful  numerical  solution  of  the  Navier-Stokes  equations  can 
yield  an  accurate  prediction  of  the  decay  and  movement  of  a  pair  of  two-dimensional  unsteady 
vortices.  They  demonstrated  the  accuracy  of  their  results  through  comparison  with  previous 
analytic  solutions.  These  results  provide  an  excellent  foundation  for  further  efforts  at  method 
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and  code  validation.  One  of  the  goals  of  the  present  work  is  to  reproduce  their  results  in  a 
three-dimensional  steady  flow  solution.  This  makes  sense  only  within  the  context  of  the 
unsteady  crossflow  analogy.  This  analogy  states  that  two-dimensional  unsteady  solutions  can  be 
good  approximations  for  three-dimensional  steady  flows  as  long  as  the  crossflow  velocities  are 
dominated  by  the  streamwise  velocity,  where  the  streamwise  velocity  is  with  respect  to  the 
reference  frame  in  which  the  flow  is  steady. 

If  a  method  or  code  can  accurately  track  a  vortex  pair,  it  should  be  able  to  handle  vortices 
in  more  general  situations.  The  code  chosen  is  DTNS3D  because  of  its  availability,  its  ease  of 
use,  and  its  successes  in  previous  applications.  The  approach  is  to  separate  the  question  of 
fundamental  vonex  tracking  capabilities  from  the  difficulties  of  turbulent  flow  modeling. 
Therefore,  the  code's  turbulent  flow  models  are  initially  turned  off  to  provide  laminar  (or 
constant  eddy  viscosity  turbulent)  flow  solutions.  The  first  computations  deal  with  an  isolated 
vortex  with  and  without  diffusion.  Comparisons  are  made  with  an  analytic  solution.  Then  the 
code  is  applied  to  the  vortex  pair,  and  comparisons  are  made  to  the  results  of  Ohring  and  Lugt.‘ 
The  computations  are  quite  straightforward  with  the  unsteady  crossflow  analogy  indicating  that 
the  three-dimensional  steady  solutions  should  mimic  the  two-dimensional  unsteady  solutions 
quite  closely.  However,  in  three  dimensions  there  is  a  third  momentum  equation,  and  hence  a 
spatial  variation  in  the  streamwise  velocity  which  must  be  handled  by  the  flow  solver.  Thus  a 
byproduct  of  this  study  is  information  gained  about  streamwise  velocity  distributions  in  vortices 
and  the  effect  of  such  distributions  on  the  validity  of  the  unsteady  crossflow  analogy. 


FLOW  CODE  DESCRIPTION 

The  DTNS  codes  were  originally  developed  by  Joseph  Gorskl  at  NSWCCD  and  have  more 
recently  been  further  developed  and  improved  by  Gorski  and  others.  Details  of  the  mathematical 
formulation  on  which  they  are  based  and  of  the  numerical  methods  used  in  the  codes  have  been 
presented  numerous  times  elsewhere.^'^  Therefore,  only  a  brief  overview  of  the  methods  and 
features  of  the  codes  are  presented  here.  The  Navier-Stokes  equations  for  steady  incompressible 
flow  are  solved  using  pseudo-compressibility  where  an  artificial  time  dependence  is  added  to  the 
continuity  equation.  The  equations  are  then  marched  in  time  until  a  steady  state  solution  is 
achieved.  When  the  steady  state  is  reached,  all  time  derivative  terms  vanish  and  a  solution  of 
the  equations  for  steady  incompressible  flow  is  obtained.  The  equations  are  discretized  with  the 
finite-volume  formulation.  Central  differencing  is  applied  to  the  viscous  terms,  but  the  high 
accuracy  upwind  total  variation  diminishing  (TVD)  scheme  developed  by  Chakravarthy  and 
Osher^  is  used  on  the  convective  terms.  The  Jacobian  matrix  of  fluxes  is  used  to  create  a  set  of 
eigenvectors  and  eigenvalues  where  fluxes  are  computed  at  the  centers  of  cell  surfaces.  The 
upwind  scheme  divides  the  flux  between  neighboring  cells  into  parts,  each  associated  with  an 
eigenvalue.  Based  on  the  signs  of  the  eigenvalues,  the  fluxes  are  differenced  in  the  appropriate 
direction  to  achieve  unconditional  stability.  A  third-order  accurate  scheme  is  chosen.  Flux 
limiting  is  available  to  prevent  oscillations  from  creeping  into  the  solution  during  startup,  but 
flux  limiting  is  generally  not  needed  for  these  incompressible  flow  computations.  The  codes 
include  both  Baldwin-Lomax  and  k-e  turbulence  models.  When  the  k-e  model  is  used,  the  TVD 
upwind  differencing  is  also  used  on  the  k-e  equations. *  For  laminar  flow  computations  the 
turbulence  models  can  be  turned  off.  All  viscous  terms  can  be  excluded  if  solution  of  the  Euler 
equations  is  desired. 

The  equations  resulting  from  the  upwind  differencing  are  solved  in  an  implicit  coupled 
manner  with  the  use  of  either  approximate  factorization  or  an  LU-ADI  solver.  The  implicit  pan 
of  the  equations  is  discretized  with  a  first-order  accurate  upwind  scheme  for  the  convective 
terms.  This  creates  a  diagonally  dominant  system  that  requires  the  inversion  of  block  tridiagonal 
matrices.  Since  the  implicit  terms  are  zero  upon  convergence,  the  solution  retains  the  high-order 
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accuracy  of  the  explicit  part  of  the  equations.  A  spatially  varying  time  step,  based  on  the  local 
eigenvalues,  is  employed. 


MATHEMATICAL  FORMULATION 

Boundary  value  problems  are  formulated  with  a  rectangular  {x,y,z)  coordinate  system  as 
required  by  the  flow  code.  For  the  vortex  pairs,  all  lengths  are  referred  to  their  initial 
separation  ag.  For  an  isolated  vortex  this  length  scale  is  not  present,  and  a  different  length 
would  typically  be  used  for  nondimensionalization.  But,  in  order  to  facilitate  comparison  of  the 
results  for  the  isolated  vortices  with  those  for  the  vortex  pairs,  Og  is  used  for  both  cases.  The 
characteristic  velocity  in  both  cases  is  where  xj,  is  the  vortex  strength.  Thus,  in 
nondimensional  form,  the  vortices  will  initially  have  a  tangential  velocity  of  one  at  a  radial 
distance  of  one  from  their  center,  and  the  initial  z  (upward)  velocity  of  the  vortex  pairs  will  be  one. 
Calculations  are  carried  out  in  parallelopiped  boxes  extending,  in  the  streamwise  direction,  from 
an  upstream  face  at  x  =  0  to  a  downstream  face  at  x  =  {  where  i  can  be  varied  depending  on  how 
long  a  trajectory  segment  it  is  desired  to  compute.  The  vertical  box  sides  are  at  y  =  0  and  y  =  20  for 
the  vortex  pair  calculations  but  the  box  contains  -20  <  y  <  20  for  calculations  with  a  single  vortex. 
For  the  pair  computation,  the  bottom  of  the  box  slants  upward  from  z  =  -10  at  x  =  0  to  z  =  -10 
+(0. 1)(  at  X  =  /.  Similarly,  the  top  slants  upward  with  the  same  slope  as  the  bottom  from  z  =  20  at 
X  =  0.  This  slant  is  included  to  accommodate  the  upward  translation  with  downstream  distance  of 
the  vortex  pair.  For  the  single  vortex  the  top  and  bottom  of  the  computational  box  are  horizontal. 

A  larger  space  is  included  in  the  box  above  the  vortices,  which  are  centered  at  z  =  0  for  x  =  0, 
because  of  some  concern  that  the  flow  at  the  slanted  top  boundary  has  a  normal  component 
directed  into  the  flow  region  and  hence  this  boundary  may  behave  as  an  inflow  boundary.  No 
problems  were  noted  with  this  boundary  for  any  of  the  computations  and  it  is  believed  that  this 
precaution  is  unnecessary. 

Within  this  box  region,  the  Euler  or  Navier-Stokes  equations  are  solved  depending  on 
whether  it  is  desired  to  consider  convective  effects  alone  or  to  include  diffusive  effects  as  well. 

This  solution  is  accomplished  through  execution  of  the  DTNS3D  code  with  turbulence  models 
turned  off  to  get  a  Navier-Stokes  solution  or  with  all  viscous  terms  turned  off  to  get  an  Euler 
solution.  The  dependent  variables  used  by  this  code  are  u,  v,  w,  and  p,  the  velocities  in  the  x,  y,  and 
z  directions  and  the  pressure,  respectively.  For  this  study  these  variables  are  interpreted  as 
dimensionless  quantities.  At  the  inflow  boundary,  the  three  velocity  components  are  specified 

(w,v,w)  =  (U(y^),V(y^),W(y^))  at  x  =  0  (1) 

At  the  plane  of  symmetry,  for  the  case  of  a  vortex  pair  the  symmetry  conditions 

—  =  0,  —  =  0  and  v  =  0  af  y  =  0  (2) 

dy  dy 

are  applied.  At  all  other  boundaries,  the  characteristic  directions  are  such  that  simple  first-order 
extrapolation  can  be  applied  to  the  velocities  as  a  numerical  approximation  to: 


d(u,v,w)  ^  Q 
dn 


(3) 
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where  n  represents  the  coordinate  direction  normal  to  the  boundary.  Pressure  boundary  conditions 
are  not  needed  for  the  mathematical  problem.  For  computational  purposes  first  order  extrapolation 
is  applied  to  the  pressure  at  all  boundaries: 


^  =  0. 

dn 


(4) 


The  pressure  is  fixed  at  p  =  1  at  one  upstream  comer  of  the  box  to  eliminate  the  arbitrariness  in  the 
pressure  definition. 

Because  this  is  a  steady-state  formulation,  initial  conditions  are  not  needed  for  the 
mathematical  statement  of  the  problem.  However,  the  numerical  approach  is  iterative,  and  initial 
conditions  are  needed  for  the  computations.  Two  basic  initial  conditions  are  utilized.  In  the  first, 
the  flow  field  is  of  constant  pressure  everywhere  with  uniform  stream  throughout  except  at  the 
upstream  boundary  where  the  desired  inflow  velocity  distribution  is  specified.  Vortices  trailing  a 
body  will  be  steady  in  a  frame  of  reference  translating  with  the  body/vortex  system.  The  uniform 
stream  referred  to  is  this  translational  velocity  and  is  the  velocity  in  the  far  field  relative  to  this 
reference  frame.  In  the  second  form  of  the  initial  conditions,  a  distribution  of  velocities  and 
pressure  is  specified  throughout  the  computational  region.  This  distribution  is  selected  to  try  to 
better  approximate  the  anticipated  steady  solution.  The  first  form  of  initial  conditions  leads  to  large 
initial  streamwise  gradients  at  the  inflow  boundary  followed  by  a  gradual  spread  of  this 
disturbance  throughout  the  region  until  a  steady-state  solution  is  achieved.  The  second  form  of 
initial  conditions  leads  to  a  less  dramatic  adjustment  of  the  initial  vortex  to  its  steady-state 
position,  the  details  depending  on  the  quality  of  the  initial  distribution  imposed.  The  latter  has  the 
potential  of  leading  to  a  numerical  solution  more  rapidly  while  the  former  is  more  representative  of 
the  nature  of  the  initial  conditions  encountered  in  practical  applications.  In  preliminary 
computations,  it  is  verified  that  the  two  initial  conditions  lead  to  essentially  identical  steady  flow 
solutions. 


GRIDDBMG 

The  vortices  are  close  to,  or  exactly  axisymmetric  in  the  crossflow  plane;  therefore,  the  most 
efficient  use  of  grid  points  would  result  from  a  grid  based  on  polar-type  coordinates  about  the 
vortex  center.  However,  in  applied  calculations  it  is  often  necessary  to  track  vortices  that  move  on 
an  unknown  trajectory  through  a  grid  that  is  tailored  to  a  body  or  other  geometrical  feature.  For 
this  work  a  grid  is  used  that  is  rectangular  in  the  crossflow  plane  to  avoid  excluding  some  of  the 
difficulties  that  might  be  encountered  through  use  of  a  grid  which  is  not  conformed  to  the  vortex 
shape.  However,  a  refinement  is  built  into  this  grid  to  assure  that  increased  resolution  is  present  in 
the  vortex  cores  at  least  for  some  distance  downstream  from  the  inflow  boundary.  A  typical  portion 
of  a  grid  cross-sectional  plane  for  the  vortex  pair  is  displayed  in  Figure  1.  The  minimum  spacings 
are  6y  =  0.0033  and  6z  =  0.0025  centered  at  y  =  0.5  and  z  =  0.0.  The  grid  spacing  expands 
logarithmically  in  both  y  and  z  directions.  For  an  isolated  vortex  the  grid  extends  to  y  <  0  and  is 
symmetric  about  y  =  0.5.  A  side  view  of  an  x,z  grid  plane  at  y  =  0.5  for  a  vortex  pair  computation 
is  shown  in  Figure  2.  This  grid  slopes  upward  parallel  to  the  top  and  bottom  boundaries.  Thus  the 
region  of  fine  resolution  tends  to  follow  the  upward  slope  of  the  vortices.  However,  to  the  extent 
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Fig.  1.  Typical  grid  in  the  y,z  plane. 


X 

Fig.  2.  Grid  cross  section  in  the  x,z  plane  at  y  =  0.5. 
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that  the  upward  slope  of  the  vortices  diminishes  due  to  diffusion  and  three-dimensional  effects, 
they  will  tend  to  drift  away  from  the  minimum  spacing  with  distance  x  from  the  inflow  boundary. 
The  basic  spacing  in  the  streamwise  direction  is  =  0.25.  For  the  single  vortex  there  are  21 
planes  separated  by  this  spacing,  extending  from  jc  =  0.0  to  x  =  5.0,  which  comprise  the  main 
region  of  interest.  These  are  followed  by  another  10  planes  where  each  succeeding  spacing  in  the 
x:-direction  is  increased  by  10  percent  over  the  previous  spacing.  These  latter  planes  are  used  to 
remove  the  outflow  boundary  from  the  region  of  interest  and  thus  minimize  any  deleterious  effects 
of  outflow  boundary  conditions.  This  box  for  the  isolated  vortex  extends  from  x  =  0.0  almost  to 
=  9.0  and  is  referred  to  as  the  short  box.  For  vortex  pair  computations,  the  region  of  interest 
extends  for  101  planes  to  jc  =  25.0.  This  is  again  followed  by  a  10-plane  region  of  expanding 
x-spacing,  which  extends  the  box  almost  to  x  =  29.0.  This  box  is  referred  to  as  the  long  box.  In 
addition  the  vortex  pair  calculations  are  continued  successfully  through  a  similar  box,  with  a 
region  of  interest  from  x  =  25.0  to  =  50.0,  by  using  the  computed  result  a.tx  =  25.0  as  an  inflow 
condition. 

There  are  a  total  of  31x81  x1 49  =  374,139  points  in  the  x-,  y-,  and  z-directions, 
respectively,  for  the  single  vortex  computations  and  1 1 1  x  61  x  149  =  1,008,879  points  for  the 
vortex  pair  computations. 


RESULTS 

SINGLE  VORTEX  WITH  NO  DIFFUSION 

An  isolated,  axisymmetric  vortex  in  the  absence  of  diffusion  should  continue  indefinitely  in 
the  streamwise  direction  with  unchanged  crossflow  characteristics.  In  order  to  test  the  ability  of  the 
code  to  handle  this  rather  simple  case,  the  viscous  terms  are  turned  off  and  the  code  is  run  in  its 
Euler  equation  mode.  The  crossflow  velocities  imposed  at  the  inflow  boundary  correspond  to 
Lamb's  two-dimensional  unsteady  (decaying)  vortex  of  initial  strength  x-where  the  vorticity  o)  is 
given  by 


6)  =  iRe/2t)e-‘^"‘^‘, 


(5) 


evaluated  at  time  t  =  2  and  for  Reynolds  number  Re  =  x/v  =  100;  where  v  is  the  kinematic 
viscosity,  and  =  (y-0.5)^  -i-  z^.  Because  the  desired  solution  has  all  A:-derivatives  zero,  it 
essentially  should  be  independent  of  the  streamwise  velocity.  A  streamwise  velocity  of  m  =  =  10 

is  imposed  at  the  inflow  plane  in  anticipation  of  the  more  complicated  cases  considered  below. 

The  code  is  run  in  the  small  box  with  the  above  inflow  condition  and  with  two  initial 
conditions,  each  for  somewhat  over  1000  iterations.  Each  iteration  consumed  about  four  minutes 
on  a  Convex  3840  computer.  For  an  initial  condition  consisting  of  the  exact  solution,  namely  the 
dependent  variables  corresponding  to  Equation  (5)  extended  throughout  the  x:-direction, 
longitudinal  cross-sectional  plots  of  the  vorticity  contours  at  the  vortex  center  are  presented  in 
Figure  3.  Even  with  the  highly  expanded  vertical  scale  relative  to  the  horizontal  scale,  the 
computed  contours  agree  quite  well  with  the  exact  (horizontal)  contours.  This  result  shows  that  the 
code  is  capable  of  producing  an  excellent  approximation  to  the  exact  solution  for  this  simple  case, 
at  least  when  the  initial  condition  is  very  close  to  the  desired  result.  To  better  test  the  ability  of  the 
code  to  generate  a  solution  from  scratch,  this  case  was  rerun  with  a  uniform  stream  (zero 
crossflow)  initial  condition.  The  extra  work  needed  to  develop  the  vortex  throughout  the  box  is 
reflected  by  the  fact  that  the  results  for  the  first  initial  condition  were  essentially  unchanged  from 
iteration  200  to  iteration  1000,  while  the  computations  for  this  second  initial  condition  took  about 
three  times  as  long  to  converge.  After  about  600  iterations  the  results  are  essentially  the  same  as 
those  presented  in  Figure  3. 
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Fig.  3.  Computed  and  exact  axial  vorticity  contours 
at  y  =  0.5  for  a  single  vortex  with  no  diffusion. 

SINGLE  VISCOUS  VORTEX  WITH  UNIFORM  AXIAL  INFLOW 

An  isolated,  axisymmetric  vortex  in  the  wake  of  flow  past  a  body  will,  in  reality,  decay  with 
distance  from  the  body  because  of  the  effects  of  molecular  or  turbulent  diffusion.  This  process  can 
be  handled  in  the  laminar  flow  case  by  the  Navier-Stokes  equations,  which  can  also  be  appropriate 
for  turbulent  flow  if  a  constant  eddy  viscosity  model  is  accurate.  To  test  the  capability  of  the  code 
to  handle  such  a  decaying  vortex,  the  computation  of  the  above  section  is  repeated  with  the  viscous 
terms  turned  on  so  as  to  run  the  code  in  its  Navier-Stokes  mode.  For  this  steady  three-dimensional 
problem,  there  is  no  exact  solution  as  in  the  previous  case.  However,  Lamb's  exact  solution  to  the 
unsteady  decay  of  a  two-dimensional  vortex  provides  a  useful  approximate  solution  as  long  as  the 
unsteady  crossflow  analogy  is  valid.  This  analogy  holds  when  the  streamwise  (axial)  velocity  is 
suitably  large  compared  to  the  crossflow  velocity  components.  This  condition  should  be  fulfilled 
as  long  as  the  axial  velocity  does  not  drop  much  below  its  inflow  value  {Uq  =  10).  Then  the 
crossflow  fields  become  essentially  independent  of  the  axial  flow  field,  and  the  three-dimensional 
steady  solution  can  be  usefully  approximated  by  a  corresponding  two-dimensional  unsteady 
solution.  The  inflow  condition  for  the  desired  three-dimensional  problem  must  be  used  as  the 
initial  condition  for  the  two-dimensional  problem.  Since  the  inflow  condition  used  for  the  current 
study  is  Lamb’s  solution  att  =  to=2,  the  complete  two-dimensional  unsteady  solution  will  be 
given  by  Equation  (5).  The  approximate  three-dimensional  steady  solution  is  then  obtained  by 
replacing  t  in  Lamb's  solution  (Equation  5)  according  to 

(6) 
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The  approximate  solution  obtained  in  the  above  manner  is  used  as  the  initial  condition  for  the 
computation.  After  a  few  hundred  iterations  a  converged  solution  is  obtained.  The  computed 
solution  must  agree  with  the  analytic  solution  near  the  inflow  plane  because  of  the  inflow 
boundary  conditions.  However,  the  computed  solution  will  agree  with  the  analytic  solution  farther 
downstream  only  to  the  extent  that  the  code  accurately  reproduces  the  diffusion  process  of  the 
governing  equations.  To  assess  the  accuracy  of  the  code,  comparisons  of  the  two  solutions  are 
made  in  Figures  4  and  5.  Figure  4  displays  a  comparison  of  axial  vorticity  contours  in  a 
longitudinal  plane,  and  Figure  5  displays  pressure  contours  in  a  cross  section  at  x  =  5.  These  two 
comparisons  display  excellent  agreement  of  the  computed  solution  with  the  approximate  analytic 
solution.  This  good  agreement  demonstrates  both  the  accuracy  of  the  code  and  the  validity  of  the 
unsteady  crossflow  analogy. 


Fig.  4.  Computed  and  analytic  axial  vorticity  contours 
at  y  =  0.5  for  a  single  vortex. 
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Fig.  5.  Computed  and  analytic  pressure  contours 
at  X  =  0.5  for  a  single  vortex. 


The  inflow  axial  velocity,  =  10  ,  was  chosen  to  assure  the  validity  of  the  unsteady 
crossflow  analogy.  It  is  possible  to  explore  the  limitations  of  this  validity  by  imposing  different 
inflow  values.  To  this  end  the  computations  are  repeated  with  an  inflow  condition  of  Uo  =  5. 

Figure  6  shows  the  resulting  axial  vorticity  contours  in  the  vortex  centerplane  compared  with  the 
analytic  solution  based  on  the  unsteady  crossflow  analogy.  Note  that  the  length  of  the 
computational  box  has  been  halved  in  correspondence  with  the  halving  of  the  axial  velocity.  It  is 
seen  that  the  close  agreement  of  Figure  4  has  been  lost.  Evidently  three-dimensional  effects  are 
large  enough  that  the  unsteady  crossflow  analogy  is  less  accurate.  It  appears  that  the  vortex 
dissipates  more  rapidly  than  for  the  case  with  larger  axial  velocity.  Actually  the  dissipation  is  the 
same,  but  axial  vorticity  is  converted  to  tangential  vorticity.  The  generation  of  tangential  vorticity 
is  demonstrated  in  Figure  7,  a  plot  of  y-vorticity,  which  corresponds  to  tangential  vorticity  in  the 
vertical  vortex  centerplane.  These  three-dimensional  steady  results  show  tangential  vorticity 
increasing  in  the  jc-direction  throughout  the  computational  box.  The  unsteady  crossflow  solution 
does  not  account  for  such  conversion  of  vorticity  from  one  component  to  another. 

Figure  8  displays  a  comparison  of  pressure  contours  in  the  cross  section  at  x  =  2.5.  This 
comparison  can  be  contrasted  directly  with  that  of  Figure  5  for  the  larger  axial  flow.  In  spite  of  the 
deterioration  with  decreasing  axial  flow  of  the  unsteady  crossflow  analogy,  it  still  provides  what 
might  be  a  useful  approximation  for  many  purposes.  It  is  not  possible  to  draw  a  general  conclusion 
on  a  precise  limit  on  the  utility  of  this  approximation. 
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Fig.  6.  Computed  and  analytic  axial  vorticity  contours 
at  y  =  0.5  for  a  vortex  with  L/j,  =  5. 


Fig.  7.  Computed  y-vorticity  contours 
at  y  =  0.5  for  a  vortex  with  Ug  =  5. 


>’ 

Fig.  8,  Computed  and  analytic  pressure  contours 
at  X  =  2.5  for  a  vortex  with  Ug  =  5. 


SINGLE  VORTEX  WITH  NONUNIFORM  AXIAL  INFLOW 

Even  when  the  unsteady  crossflow  analogy  is  valid,  there  is  a  third  momentum  equation  (for 
axial  momentum)  that  must  be  satisfied  in  the  three-dimensional  case.  Since  the  pressure 
distribution  changes  as  the  vortex  decays,  there  are  streamwise  pressure  gradients.  These  appear  in 
the  axial  momentum  equation,  and  hence  they  must  be  balanced  by  gradients  of  the  axial  velocity. 
These  gradients  affect  the  crossflow  solution  through  the  other  momentum  equations  but  only 
negligibly  as  long  as  the  conditions  of  the  crossflow  analogy  are  not  violated.  The  computed  axial 
velocity  contours  at  x  =  5  for  the  flow  of  the  previous  section  with  inflow  Ug  =  10  are  displayed  in 
Figure  9.  The  deviation  from  the  inflow  velocity  is  around  one  percent.  The  deviation  of  the 
contours  from  axisymmetric  is  due  to  numerical  error,  which,  not  surprisingly,  is  apparent  at  these 
low  axial  velocity  variation  levels. 

It  would  seem  that  any  inflow  distribution  of  axial  velocity  can  be  specified  and  a  suitable 
Navier-Stokes  solution  can  be  generated.  If  the  axial  velocity  dominates  the  crossflow  velocities 

solution  and  will  be  independent  of  the  axial  inflow.  Batchelor’  derived  an  approximate  solution 
for  the  axial  flow  distribution  based  on  the  crossflow  analogy.  According  to  this  solution  the  axial 
flow  is  given  by 

u  =  C/„-(Re/&c)log(C4Rex)e,(ti).(Re/8x)e/ti)-C(Rc/8x)e,(q),  (7) 

where 

Tj  =  ReU/V^x,  (8) 
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Equation  7  is  used  to  verify  the  code's  ability  to  accurately  compute  the  axial  velocity 
distribution.  For  Uq  =  10  that  equation  yields  m  =  [/(,=  10  far  from  the  vortex  core.  The  constant  C 
is  chosen  to  yield  a  desired  value  of  u  at  the  center  of  the  vortex,  r  =  0,  at  the  inflow  plane,  x  =  0. 
Initially,  C  is  chosen  such  that  the  inflow  axial  velocity  falls  10  percent,  from  its  farfield  value  to  9 
at  the  vortex  center.  The  resulting  inflow  axial  velocity  distribution  is  shown  in  Figure  10. 
Corresponding  velocity  contours  at  x  =  5  are  shown  in  Figure  1 1 .  It  is  seen  that  both  the 
computations  and  the  approximate  analytic  solution  display  an  increase  in  the  streamwise  direction 
of  the  axial  velocity  at  the  center.  The  analytic  solution  shows  an  increase  from  9.0  at  x  =  0  to 
almost  9.1  at  X  =  5,  while  the  numerical  results  show  an  increase  to  a  value  somewhat  higher  than 
9.1.  Apart  from  such  slight  differences,  the  computed  and  analytic  solution  agree  quite  well.  The 
computed  and  analytic  vorticity  and  pressure  contours  agree  even  more  closely  than  the  axial 
velocities  with  no  difference  noticeable  between  the  two  solutions  when  they  are  plotted  on  the 
scale  of  Figure  11.  The  small  differences  in  the  axial  velocity  contours  are  attributed  to  a  slight 
breakdown  of  the  unsteady  crossflow  analogy  due  to  the  axial  velocity  variations,  and  hence  show 
the  limitations  of  the  approximate  analytic  solution.  This  explanation  is  verified  by  repeating  the 
computations  for  an  inflow  axial  velocity  which  varied  only  1  percent  from  the  farfield  to  the 
vortex  center,  in  which  case  the  axial  velocity  contours  from  the  two  solutions  at  x  =  5  fall  on  top 
of  each  other  when  plotted  as  in  Figure  11. 


Fig.  10.  Axial  velocity  distribution  imposed 
at  X  =  0  for  a  single  vortex. 
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Fig.  11.  Computed  and  analytic  axial  velocity  contours  at  x  =  5 
for  a  vortex  with  inflow  corresponding  to  Fig.  10. 

It  is  possible  to  explore  the  effects  of  large  axial  flow  variations  through  the  computational 
approach  beyond  the  limits  of  usefulness  of  the  analytic  solution.  To  this  end  computations  are 
carried  out  with  an  inflow  axial  velocity  distribution  that  varies  100  percent  from  10  in  the  farfield 
to  0  at  the  center.  This  inflow  distribution  is  shown  in  Figure  12.  Computed  and  analytic  velocity 
contours  at  x  =  5  are  plotted  in  Figure  13.  The  difference  between  the  solutions  is  much  larger  here 
than  for  the  previous  case  of  smaller  axial  flow  gradient.  The  axial  velocity  minimum  at  the  center 
is  about  1 .9  according  to  the  analytic  solution  but  about  4.4  in  the  numerical  results.  The  crossflow 
field  is  also  different  at  x  =  5  for  the  two  solutions.  Figure  14  shows  a  longitudinal  section  of  axial 
vorticity  contours  for  the  computations  and  the  analytic  solution.  The  two  flow  fields  agree  atx  = 

0,  and  both  solutions  start  there  with  a  maximum  vorticity  of  about  25.  However,  the  maximum 
grows  quickly  with  downstream  distance  x  to  a  value  of  almost  42  in  the  computations,  while  the 
analytic  solution  does  not  exhibit  this  growth.  Plots  of  the  y  component  of  vorticity  are  shown  in 
Figure  15.  This  component  is  rather  large  at  the  inflow  due  to  the  axial  velocity  gradient  there.  The 
computed  results  show  a  rapid  loss  of  y- vorticity  with  downstream  distance.  Evidently  tangential 
vorticity  is  rapidly  converted  to  axial  vorticity  for  a  short  distance  downstream  from  the  inflow 
plane.  Clearly,  this  case  is  well  beyond  the  range  of  usefulness  of  the  unsteady  crossflow  analogy 
which  accounts  for  convection  and  diffusion  of  the  two  vorticity  components  but,  as  mentioned 
previously,  cannot  account  for  transfer  between  the  components. 
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Fig.  12.  Axial  velocity  distribution  imposed 
at  X  =  0  for  a  single  vortex. 


Fig.  13.  Computed  and  analytic  axial  velocity  contours  at  x  =  5 
for  a  vortex  with  inflow  corresponding  to  Fig.  12. 


VORTEX  PAIR  WITH  UNIFORM  AXIAL  INFLOW 

A  vortex  pair  will  exhibit  the  convective  effect  of  each  vortex  on  the  other  in  addition  to  the 
effects  of  diffusion.  To  test  the  ability  of  the  code  to  cope  with  both  effects  simultaneously  the 
code  is  run  in  its  laminar  flow  mode  with  Re  =  100  and  with  the  Lamb  vortex  pair  at  =  2  and 
uniform  axial  velocity  as  the  inflow  condition.  Calculations  have  been  carried  out  in  the  long  box. 
After  convergence  is  achieved,  the  results  at  x  =  25  are  used  as  inflow  conditions  to  a  similar  box 
to  continue  the  solution  to  about  x  =  50  in  the  streamwise  direction.  The  continuation  of  the 
solution  in  this  manner  works  quite  well.  There  is  not  an  accurate  analytic  solution  for  comparison 
in  this  case.  However,  the  excellent  two-dimensional  unsteady  computational  results  of  Ohring  and 
Lugt'  can  be  used. 

A  plot  of  the  axial  vorticity  contours  in  the  longitudinal  centerplane  of  one  of  the  vortices  is 
presented  in  Figure  16.  The  vortex  moves  upward  because  of  convection.  At  the  same  time 
diffusion  acts  to  spread  the  vorticity  and  reduce  its  maximum.  This  results  in  a  downward 
curvature  of  the  vortex  which  increases  with  downstream  distance  and  is  most  visible  in  the  upper 
5  level  contour.  This  behavior  represents  a  reduction  in  rise-velocity  and  is  clearly  visible  in  the 
plot  in  Figure  17  where  the  approximate  vertical  location  of  the  pair  is  plotted  versus  time.  This 
rise  history  is  compared  with  the  results  of  Ohring  and  Lugt  and  with  the  rise  history  for  a  pair  of 
(inviscid)  point  vortices.  The  oscillations  in  the  current  results,  most  visible  at  the  later  times,  are  a 
result  of  the  rather  crude  method  for  approximating  the  vortex  location.  The  vorticity  is  computed 
at  grid  cell  centroids.  The  location  of  the  vortex  is  taken  to  be  at  the  grid  cell  centroid  with  the 
largest  vorticity  value.  The  vortex  moves  upward  with  distance  x,  as  does  the  grid.  Therefore,  the 
maximum  vorticity  tends  to  appear  in  the  same  cell  for  several  successive  constant  x  planes.  When 
the  vortex  has  moved  enough  for  the  maximum  vorticity  to  appear  in  a  different  grid  cell  at  one  x 
plane  than  in  the  previous  plane,  there  appears  a  kink  in  the  rise  history.  The  oscillations  grow 


Fig.  16.  Computed  axial  vorticity  contours  at  y  =  0.5 
for  a  vortex  pair. 
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Fig.  17.  Vertical  location  of  vortex  pair  versus  time. 


with  time  since  the  vortices  gradually  move  away  from  the  area  of  maximum  grid  resolution, 
which  was  chosen  to  follow  the  point  vortex  pair  trajectory.  A  smoother  rise  history  could  be 
generated,  if  warranted,  by  better  estimating  the  location  of  the  vorticity  maximum  within  a  grid 
cell.  The  current  results  predict  a  rise-velocity  that  is  slightly  lower  than  that  of  Ohring  and  Lugt. 
Possible  explanations  for  this  include  numerical  inaccuracies,  differences  in  starting  conditions, 
and  real  differences  between  the  two-dimensional  and  three-dimensional  flows. 

To  assess  the  impact  of  grid  spacing  on  the  results,  the  first  long  box  computation  is  repeated 
with  every  other  grid  point  removed  in  all  three  directions.  Except  for  a  further  loss  of  smoothness 
the  resulting  rise-history  curve  is  essentially  unchanged,  as  seen  in  Figure  18  where  it  is  compared 
with  that  for  the  original  grid. 

To  assess  the  effect  of  starting  conditions  on  the  results,  the  first  long  box  computation  with 
the  vortex  pair  is  repeated  with  an  inflow  condition  corresponding  to  a  pair  of  Lamb  vortices  at  tg  = 
1.5.  The  resulting  rise  history  is  presented  in  Figure  19,  where  it  is  compared  with  the  previous 
results  for  to  =  2  and  with  the  results  of  Ohring  and  Lugt  for  which  tg  =  0.2.  The  current  results  for 
the  two  starting  conditions  are  quite  close.  Moreover,  the  slight  dependence  on  tg  is  in  the  wrong 
direction  to  provide  an  explanation  of  the  deviation  from  the  results  of  Ohring  and  Lugt. 
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A  third  possible  explanation  for  the  deviation  from  the  two-dimensional  unsteady  results  lies 
in  legitimate  three-dimensional  effects.  One  such  effect  is  the  generation  of  streamwise  velocity 
gradients  through  the  axial  momentum  balance.  For  example,  Figures  20  and  21  show  computed 
axial  velocity  contours  at  two  x  locations.  At  a;  =  12.5  the  computed  axial  velocity  shows  a 
minimum  between  9.6  and  9.7  near  the  vortex  centers  and  a  region  of  increased  velocity  with  a 
maximum  between  10.1  and  10.2.  By  x  =  25  the  contours  are  similar  but  the  local  maximum  has 
decreased  to  below  10.1 .  This  three-dimensional  effect  is  certainly  small  but  could  have  some 
influence  on  the  rise  velocity.  To  further  assess  this  influence,  additional  computations  are  carried 
out  with  uniform  inflow  axial  velocities  of  5  and  100.  As  already  shown  for  the  single  vortex, 
three-dimensional  effects  are  enhanced  by  a  decrease  of  axial  velocity.  Figure  22  shows  a 
comparison  of  the  computed  rise  histories  for  the  three  axial  flow  rates.  The  rise  of  the  vortices  is 
lowest  for  the  smallest  axial  inflow  (  =  5).  In  that  case,  more  axial  vorticity  is  converted  to 

tangential  vorticity  resulting  in  a  reduced  tangential  velocity  and  thus  a  reduced  upward  velocity  of 
the  vortex  pair.  On  the  other  hand,  the  largest  axial  flow  case  ( t/o  =  100 ),  exhibits  the  highest 
rise.  This  result  clearly  demonstrates  the  nature  of  the  two-dimensional  unsteady  flow  as  the  limit 
of  the  three-dimensional  steady  flow  with  infinite  axial  velocity.  It  also  shows  that  for  trajectories 
of  the  length  of  those  considered  here,  the  crossflow  analogy  yields  a  very  good  estimate  of  the 
vortex  pair  position  as  long  as  Ug  >  10. 

VORTEX  PAIR  WITH  NONUNIFORM  AXIAL  INFLOW 

To  investigate  the  effect  of  axial  velocity  variations  on  the  rise  of  the  pair,  computations  are 
carried  out  with  an  axial  velocity  inflow  given  by  Batchelor's  approximate  solution  (Equation  7) 
with  the  constant  C  chosen  such  that  the  axial  velocity  is  zero  at  the  inflow  vortex  centers.  As  for  a 
single  vortex  with  this  axial  flow  variation,  part  of  the  tangential  vorticity  inherent  in  this 
distribution  is  quickly  converted  axial  vorticity.  However,  as  shown  in  Figure  23,  this  increased 
axial  vorticity  (and  corresponding  tangential  velocity),  is  initially  confined  to  the  vicinity  of  the 
vortex  core  and  has  little  influence  on  the  rise-velocity  of  the  pair.  Later,  as  the  increased  axial 
vorticity  spreads,  the  rise- velocity  and  vertical  location  of  the  pair  noticeably  exceeds  that  of  the 
pair  with  uniform  axial  flow.  It  is  interesting  that  the  overall  reduction  of  uniform  axial  flow 
results  in  a  reduction  in  the  rise  of  the  pair  (Figure  22),  while  the  local  reduction  of  the  axial  flow 
in  the  vortex  cores  results  in  an  increase  in  the  rise  of  the  pair  (Figure  23). 
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CONCLUSIONS  AND  RECOMMENDATIONS 


A  general-purpose  RANS  computer  code  has  been  tested  for  its  ability  to  track  isolated 
trailing  vortices  and  vortex  pairs.  It  is  shown,  through  comparison  with  analytic  solutions  and 
previous  computations,  that  the  code  can  accurately  locate  such  vortices  at  least  for  the  case  of 
laminar  flow  and  for  cases  where  turbulent  flow  is  adequately  modeled  with  a  constant  eddy 
viscosity.  Computations  with  spatially  varying  inflow  axial  velocity  demonstrate  the  effect  of  this 
flow  component  including  its  impact  on  the  validity  of  the  unsteady  crossflow  analogy. 

Many  vortices  of  practical  interest  are  turbulent,  and  the  constant  eddy  viscosity  model  of 
turbulence  is  too  crude  for  their  accurate  treatment.  For  this  reason,  it  is  desirable  to  extend  this 
work  to  a  detailed  study  of  turbulent  trailing  vortices.  The  impact  of  currently  used  turbulence 
models  on  the  computational  results  can  be  assessed.  Other  promising  models  can  be  tested.  The 
absence  of  analytic  solutions  and  the  shortage  of  detailed  measured  data  will  be  a  problem.  A  goal 
will  be  the  selection  or  development  of  a  model  that  adequately  reflects  the  physics  of  these  flow 
structures.  Such  a  turbulence  model  implemented  in  a  RANS  code  with  an  accurate  vortex 
tracking  capability  such  as  demonstrated  in  this  work  will  provide  a  valuable  computational  design 
tool. 
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